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DEVELOPMENT OF COMPUTER PROGRAM NAS3D USING VECTOR PROCESSING 
FOR GEOMETRIC NONLINEAR ANALYSIS OF STRUCTURES 

By 

P. D. Mangalgiri 1 

Principal Investigator: Ramamurthy Prabhakaran 2 

ABSTRACT 

An algorithm for vectorized computation of stiffness matrices of an 8- 
noded isoparametric hexahedron element for geometric nonlinear analysis was 
developed. This was used in conjunction with the earlier 2-D program GAMNAS 
to develop the new program NAS3D for geometric nonlinear analysis. A con- 
ventional, modified Newton-Raphson process is used for the nonlinear analy- 
sis. New schemes for the computation of stiffness and strain energy release 
rates is presented. The organization of the program is explained and some 
results on four sample problems are given. The study of CPU times showed 
that savings by a factor of 11-13 were achieved when vectorized computation 
was used for the stiffness instead of the conventional scalar one. Finally, 
the scheme of inputting data is explained in the Appendix. 
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1. INTRODUCTION 

Many problems in the study of the behavior of composites require geo- 
metric nonlinear analyses due to large strains and displacements involved. 
Examples of these are the behavior of laminated composites under compression 
(buckling), with or without delamination, and deformation of lap shear geo- 
metries. Often, the two-dimensional (2-D) analysis is used successfully to 
study stress states and estimate fracture mechanics parameters in these 
problems. However, such an analysis proves inadequate in many cases where 
there exist complexities in the composite lay-up, geometry or form of delam- 
ination. A full three-dimensional (3-D) analysis is then required. Such an 
analysis, and particularly the nonlinear one, involves extensive computa- 
tions and is normally prohibitively time-consuming with conventional comput- 
ers and algorithms. The advent of vector processor computers, like the VPS- 
32 system at NASA Langley, has provided an economically feasible solution to 
the time-consuming operations of 3-D analyses. However, to fully utilize 
the fast processing capability of these computers in the vector form, new 
algorithms are required. These algorithms should be aimed specifically at 
taking advantage of vector processing by carrying out computations by build- 
ing long vectors. 

The two important time-consuming stages in the Finite Element (FE) 
analysis of structures are the computation of stiffness matrix and solution 
of resulting system of equations. The development of solution routines has 
received considerable attention in the past and already efficient vectorized 
solution routines are in operation on the VPS-32 system (Refs. 1, 2). Com- 
putation of stiffness in the vector form was studied by Noor and Hartley 
(Ref. 3) who developed an algorithm which builds and uses long vectors. 
Savings by a factor of 5 to 6 were obtained in CPU time when routines based 



on these algorithms were used. Stiffness calculations based on these algo- 
rithms were effectively used in linear 3-0 analysis by Raju et al . (Ref. 4). 
However, this has not been applied in the case of nonlinear analysis so far. 
Since the nonlinear stiffness calculation involves many more computations 
than the linear one, and since the nonlinear solution process requires 
repeated calculation of stiffnesses, the vectorization in this segment 
should lead to decreasing computation time. The present effort is directed 
towards this goal. 

Development of a capability for general nonlinear analysis of struc- 
tures with cracks or discontinuities is a complex, but worthwile exercise. 
The present work is a step in this exercise and is basically aimed at ana- 
lyzing laminated composites with delaminations under compression and lap 
shear geometries which involve large rotations during deformation. A good 
review of many efforts to estimate stresses and fracture mechanics param- 
eters in these problems can be found in references 5-7. Almost all of this 
work except Lof's (Ref. 7) is based on 2-D FE analysis with plane stress or 
plane strain. It is now generally recognized that cracked lap-shear geom- 
etries have large rotations and require a geometric nonlinear analysis to 
get accurate estimates of strain energy release rate in various modes of 
fracture [6]. A 2-D nonlinear analysis was developed in an earlier effort 
at NASA Langley (Ref. 6) and resulted in the program GAMNAS for the case of 
plane stress and plain strain idealization (Ref. 8). Subsequently, a quasi- 
3-D linear analysis was developed for the case of a symmetric double cracked 
lap shear specimen (Ref. 9). This proved to be a good compromise between 
the accuracy of the results and the cost of the solution for the 3-D nature 
of the problem. However, this analysis involved superpositions and thus 
could not be extended to the case of nonlinear behavior. 
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In the present work, the program GAMNAS and the vectorized algorithm 
for linear stiffness computation used in earlier 3-D linear analysis are 
taken as the starting points. Routines are developed for vectorized compu- 
tation of nonlinear stiffness using 8-noded isoparametric hexahedron ele- 
ment. These are used to develop a code for the 3-D geometric nonlinear 
analysis. The overall flow of the program follows closely that of GAMNAS. 
Some theoretical aspects which served as the basis for the development of 
the present program are presented. First, the general formulation of the 
nonlinear analysis is presented. Next, schemes of computation of stiffness 
and of extracting strain energy release rates are discussed. Then, the 
description and organization of the program NAS3D is given along with the 
study of computational times involved in various stages. Finally, the 
results of some simple problems are presented as test cases and examples. 
The details of input data required are given in the Appendix. 

2. THEORY 

A displacement based FE formulation is used in the present work. Such 
a formulation of a geometric nonlinear problem is given in Ref. 10. The 
salient features are given in this section. Further, schemes for computa- 
tion of stiffness and strain energy release rates are also discussed. 

2.1 Formulation 

(i) Strain-Displacement Relations 

A general definition of strains using Green's strain tensor is used to 
take into account the large displacements and rotations. This defines the 
strains as follows: 
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other components being obtained by suitable permutation of subscripts. It 
may be noted that if the displacements are small, the general first order 
strain approximation is obtained by neglecting the quadratic terms. 


(ii) Stress-Strain Relations 

For small strains the stress-strain relations can be written in general 
form as 


io] = [D] ({e} - (e} 0 ) + (a} 0 (2) 

in which [D] is the usual set of the elastic constants and the subscript 0 
refers to initial values. When the strains are large, a nonlinear stress- 
strain relationship for the material has to be used. This is not considered 
in the present formulation. 

(ii) Equilibrium Equations 

Equilibrium conditions have to be satisfied between internal and 
external generalized forces. If the displacements are prescribed in the 
usual manner of an FE analysis by a finite number of nodal parameters {a}, 
the necessary equilibrium conditions can be obtained by using principle of 
virtual work as when { } represents the sum of internal and external 
generalized forces and [B] is defined from the strain-displacement relations 
of Eq. (1) as 
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( 3 ) 


{*(«)} = / [B] T {a} d V - {f} = 0. 

V 

Since the strains depend nonlinearly on displacements, the matrix B is now 
dependent on {a}. Also, the quadratic nature of Eq. (1) implies that [B] is 
linear with {a}. The stresses {a} are linearly related to the strains. 
Hence, Eq. (3) results in cubic simultaneous equations. 

(iv) Iterative Solution 

The governing equations resulting from Eq. (3) are solved iteratively 
using modified Newton-Raphson process (Ref. 10). The basic Newton-Raphson 
method for the first load step may be outlined as follows: 

1. Obtain a linear solution using the linear stiffness matrix [K 0 ]: 

{ ^1 = [Ko { f) ; (4) 

2. Calculate residuals W with Eq. (3); 

3. Check for convergence. Stop if (4»> is sufficiently small; 

4. Calculate tangential stiffness matrix [K^] as defined by the 
equation 

d{*} = [K t ] d{ a} ; (5) 

5. Solve for correction to displacements: 

{4a} = - [K-j-3 -1 {A*} ; (6) 

6. Update the displacements: 



{a} = {a} + {Aa} ; 


(7) 


7. Go to step 2. 

When multiple load steps are used, only step 1 changes. After obtaining a 
converged solution for load step "i," the initial solution for the next load 
step is 


<a> 1+1 - <a,) + [K t ]-' (if), (8) 

where {Af} ^ is the load increment after the i-th step. 

In actual practice a modified Newton-Raphson process is used. In this 
process, the tangential stiffness is not calculated after every iteration 
but after a predetermined number of iterations. This saves considerable 
computation time. 


2.2 Computation of Stiffness 

A linear FE analysis requires the computation of stiffness matrix [K 0 ] 
for an element in the following form (Ref. 10): 

[K 0 ] = / [B] T [D] [B] dV (9) 

V 

where [B] is the matrix of strain-displacement relations obtained from Eq. 
(1) neglecting the nonlinear terms. In a nonlinear analysis one needs to 
evaluate a tangent stiffness matrix [Ky] as defined by Eq. (6). On taking 
differentials from Eq. (3), [Ky] can be obtained as 
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[K t ] = [K 0 ] + [K l ] + [K ff ] = [K] + [K o ] (10) 

where 

W - [K 0 ] + [K l ] = /„ [B] T [D] [B] dv (11) 

and 

[K a ] = /„ [G] T [M] [G] dv. (12) 

Here, [B], as defined in Eq. (3), incorporates the nonlinear strain-dis- 
placement relationship, [G] contains the derivatives of the shape functions, 
and [M] contains the stresses in the element. The matrix [B] is obtained 
as 


[B] = [B] + [BJ (13) 

and 

[B L ] = [A] • [G] 

where [A] contains the derivatives of displacements. Thus, 
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where N's are vectors of the shape functions for the isoparametric 8-noded 
hexahedron and subscripts denote the derivatives. On carrying out multi- 
plications in Eq. (13), it can be shown that [K^J has the following struc- 
ture: 


cy • 


[H] 

[0] 



[o; 

[H] 

to. 


[0] 

[0] 

[H] 



(17) 


where [H] is defined by 


[H] = /„ [B,] T [S] [B,] dv 


where 


[Bi] = [N x , N y N z ] 


and [S] = 


a t r 
x xy xz 


T CT 

xy y 


yz 

t t a 
xz yz z 


(18) 


(19/20) 


The integrals are evaluated by numerical integration using Gaussian quad- 
rature so that Eq. (12), for example, is evaluated as 


[K] = I [B]! [0] i [B], | J| i (21) 

Note that the right-hand sides of Eqs. (10), (12), and (19) have the same 
form, and similar logic can, therefore, be used in their evaluation. Hence, 
the algorithm given by Noor and Hartley (Ref. 3) for linear stiffness compu- 
tation is also applicable in the case of the nonlinear stiffness computa- 
tion. The method used in the present work follows very closely the algo- 
rithm of Noor and Hartley, and Ref. 3 can be seen for the details. The 
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computation of [K] proceeds as follows: 


1. Long vectors of length m.n containing values of derivatives of m 
shape functions (N , N , etc.) at n integration points and products 

X Jr 

of these witn derivatives of displacements (u N , u N , etc.) are 

xx x y 

formed. These vectors are stored in arrays BV and BLV as shown in 
Fig. 1. These now contain all the values required for matrices [B] 
and [8^] in Eqs. (16, 17). In the present program, the number of 
shape functions is m=8 and 3-point integration in each coordinate 
direction is used which gives n=27. Thus, vectors in BV and BLV 
have a length of 216. 

2. The long vectors formed in step 1 are multiplied by the vector of 
weights and |J|. 

3. Array SU is formed which contains the vectors in the product [B] T 
[D]. The j-th column vector in SU is formed as the sum given by 


{SU} . . = l (BV) . * [D] . 

J ij z U £J 

where the index 1 takes the values required by the strain-displace- 
ment matrix. A pointer matrix containing applicable values of 1 is 
formed for this purpose. The array SU has the same structure as BV 
and BLV and contains vectors of length mn. (See Fig. 2). 

4. The nodal stiffness coefficients are now evaluated. For this pur- 
pose the element stiffness matrix is partitioned as shown in Fig. 

3. The independent blocks of the element stiffness matrix that 
occupy each column are formed simultaneously. The operations 
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-derivative of i-th shape function at j-th integration point 
-derivative of displacement u at j-th integration point. 
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Structure of vector SU and vector SUM replicated from SU. 
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Storage organization of the element stiffness matrix [K 




involved in the evaluation of the stiffness coefficients of the j- 

th column blocks (i.e. [K. .], [K . . ..] - - - [K. ]) are as 

1 j j , j+i J ,m 

follows: 

1. Array SUM is formed from the array SU by replicating j times 
the elements corresponding to the j-th shape function as shown 
in Fig. 2. The length of vectors in SUM is nm(ntf-l)/2. 

2. Arrays BVJ and BLVJ are formed from BV and BLV, respectively, 
by replicating the values corresponding to j-th shape function 
(m-j+1) times. Thus, the length of vectors in BVJ and BLVJ is 
now nm(m+l ) /2 . Array BV is used for the linear case. The sum 
of appropriate columns of BV and BLV as dictated by the strain- 
displacement relations is used for the nonlinear case. 

3. The matrices (K .] (l=j to m) are formed by taking the products 

X.J 

of vectors in SUM and BVJ or BLVJ. The elements of [K..] are 

^ J 

vectors of length (m-j+l)n. Once again, as in forming SU 
above, the strain-displacement relations decide which vectors 
from SUM would multiply each vector from BVJ or BLVJ. The same 
pointer matrix can be used for this purpose. 

4. Stiffness coefficients K £ j are then obtained by summing the n 
coefficients K^- - - which occupy contiguous loca- 
tions in memory. 

Calculation of [K^J proceeds in the same manner with proper definitions 
of the [B] and [D] matrices. As seen from Eq. (18), only a smaller matrix 
[H] needs to be calculated which then can be replicated to get [K r ]. 

To make the element suitable for analyzing bending deformation use is 
made of reduced integration techniques (Ref. 10). In this technique, one or 
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more of shear strains are integrated by taking a lower order representation. 
For the present case of the isoparametric hexahedron, the lower order ap- 
proximation makes the shear strains have a constant value in the element and 
can be easily incorporated in the computation of the stiffness matrix. For 
this purpose, the elements of the strain-displacement matrix which 
correspond to the strain being reduced-integrated are changed to their 
values at the centroid prior to the calculation of the stiffness. Simi- 
larly, to evaluate [K a ] , the elements which multiply the shear stress cor- 
responding to the shear strain being reduced-integrated are changed to their 
values at the centroid. 

2.3 Computation of Strain Energy Release Rates 
The strain energy release rates are calculated using a virtual crack 
extension technique similar to that reported in Ref. 11 and used in GAMNAS 
(Ref. 6, 8). This technique uses the forces transmitted across the crack 
tip to determine the energy release rates. The technique assumes that the 
change in geometry caused by a very small crack increment does not change 
the forces and displacements near the crack tip. The use of the technique 
requires that the mesh near the crack tip be uniform and synmetric about the 
crack plane. Because of the large displacements and rotations involved, the 
orientation of the crack plane is calculated in the deformed configuration 
in terms of the direction cosines of the normal to the crack increment 
plane, and the modes of fracture are redefined for this local system. The 
strain energy release rates are calculated for these modes of fracture. In 
calculating the energy, first the forces at the relevant nodes are reduced 
to force per unit length of the crack front (traction) using a strain energy 
equivalence. This is equivalent to inverting the process of obtaining the 
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consistent load vector for a given stress distribution. The energy obtained 
by multiplying these tractions by corresponding displacements is then di- 
vided by the length of the crack increment to get the strain energy release 
rates. 


3. DESCRIPTION AND ORGANIZATION OF PROGRAM NAS3D 

The program NAS3D is basically developed for 3-D geometric nonlinear 
analysis of lap shear geometries with cracks. An example of such a geometry 
would be a cracked lap shear specimen used for study of delaminations in 
composites or debonds in adhesively bonded joints. The program is written 
in FORTRAN VERSION 2 operating on CYBER 200 and VPS-32 systems at NASA 
Langley. It relies heavily on the vector processing capabilities of these 
machines for its speed and efficiency, NAS3D uses a simple 8-noded isopara- 
metric hexahedron element. The iterative scheme for the nonlinear solution 
is a modified Newton Raphson process. It outputs the nodal displacements, 
stresses at the centroids of the elements, and reactions at the nodes. For 
the cracked configuration, it outputs the strain energy release rates along 
the crackfront in all the three modes as well as the total mean strain 
energy release rate. 

The overall flow of NAS3D is similar to that of GAMNAS (Ref. 8). Even 
the material nonlinearity features of the GAMNAS are retained in the 3-D 
program, although this part has not been tested so far. This is done with a 
view to incorporate the material nonlinearity in the future. The flowchart 
of the program is shown in Fig. 4. Only one proportional load vector is 
input. The different load numbers (LOADNUM) refer to the scale factor by 
which the load vector is multiplied. For each new load, a linear incre- 
mental solution is obtained in the main program and then iterated in the 
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IMPOSED 
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ITERATE FOR NONLINEAR 
SOLUTION 

INCREMENT PROPORTIONAL 
LOAD VECTOR 


Fig. 4, Flowchart for main program. 
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routine ITERATE to obtain the incremental nonlinear solution. A restart 
option is included to facilitate the starting of the program from a solution 
to an itermediate load step. A description of the input data required for 
the program is given in the Appendix. 

Table 1 gives the computation times obtained for various computations 
with NAS30. Comparative times using conventional scalar algorithsm are also 
given. The vector version of calculating stiffness matrix reduces the CPU 
time required by a factor 11 to 13 depending upon the choice of selective 
reduced integration. It may be noted here that the calculation of residuals 
and stresses is as time consuming as the stiffness calculation. Future 
efforts may be directed towards reducing the computations or computational 
times for the residuals. 


4. SAMPLE PROBLEMS 

In this section, some results obtainedon four sample problems using 
NAS3D are given. These were used as test cases for the program and serve as 
useful examples. 


4.1 Slender Beam 

A slender beam of width 1", thickness .2" and length 20" is analyzed. 
The beam is supported at one end (see Fig. 5) and a transverse displacement 
w is applied at the center of the other end. The mesh and the configuration 
are shown in Fig. 5a. Although, the beam has no transverse stiffness ini- 
tially (at no load), the geometric nonlinear effects stiffen the system as 
the transverse displacement increases. Figure 5b shows the calculated axial 
stress at the supported end of the rod at the mid-width and the edge. The 
results of 3-0 FE analysis are shown by symbols and an exact solution ob- 
tained by the use of simple trigonometry for a rod with no bending stiffness 
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Table 1. CPU Times for Various Stages in the Nonlinear Analysis 


Nodes = 8418, Elements = 6520, Bandwidth = 456, DOF = 25254 



Operation 


Time (seconds) 

1 . 

Data Input 


16.77 

2. 

Linear Stiffness and Assembly 




(Scalar) 

(Vector) 

71.00 


Stiffness Time/Element 

(Vector) 

.0096 

3. 

Nonlinear Tangent Stiffness and Assembly 




(Scalar) 

2180.00 



(Vector) 

171.00 


Stiffness Time/Element 

(Scalar) 

0.3200 



(Vector) 

0.0251 

4. 

Stresses and Reactions 

(Vector) 

140.0 

5. 

Residuals 

(Vector) 

115.0 

6. 

Solution of Equations 

(Scalar) 

3933.0 



(Vector) 

57.0 
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is shown by the solid line. As seen from the figure, the FE results predict 
the nonlinearity very well. The 3-D results are close to the exact at the 
mid-width and only slightly lower at the edge. Table 2 gives the comparison 
of the stresses obtained with full and reduced integration. It is observed 
that at the low loads the reduced integration indicates some error. This 
shows that the reduced integration should be used judiciously and only where 
bending is involved. 


4.2 Double Cantilever Beam (DCB) 

Double cantilever beamgeometries are extensively used in the study of 
fracture of fracture of composites and adhesive bonded joints. Figure 6 
shows a symmetric DCB specimen with the FE mesh for the upper half which was 
analyzed. A load of 40 lb uniformly distributed over the 1" width was 
applied. The linear and nonlinear solutions were found to' differ by an 
insignificant amount; for example, the error in tip deflection in the linear 
solution was less than 1%. The analysis was conducted using full integra- 
tion, reduced integration of one shear strain (IREDIN1), and reduced inte- 
gration of all shear strains (IREDU=3). The tip deflection w and average 
strain energy release rate G are compared with the strength of materials 
solution in Table 3. It is observed that the full integration performs 
poorly when bending is involved. The selective reduced integration allows 
the user to choose the proper reduced integration scheme suitable for a 
given problem. The difference in the two cases of reduced integration is 
not large here (Table 3) but it could be significant if the beam deforms in 
the width direction as well. In such a case, reduced integration of all 
shear strains may prove beneficial and necessary. 

Figure 6b indicates the distribution of the strain energy release rate 
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Table 2. Example of Slender Beam Axial Stress in the Beam at Various 
Deflections. 


Deflection 
w in 

Axi< 

il stress (psi) 


Reduced 

Integration 

Full Integration 

Theory 

Elm. 1 

Elm. 2 

Elm. 1 

Elm. 2 

1 

9.322 

9.373 

12.45 

12.41 

12.49 

2 

44.66 

44.89 

49.72 

49.54 

49.88 

3 

105.77 

106.49 

111.49 

111.06 

111.87 

5 

299.96 

303.35 

307.36 

306.02 

307.76 

7 

589.13 

596.16 

598.8 

595.9 

594.81 
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Example of double cantilever beam, (a) con 
of G along the crack front. 


Table 3. Example of Double Cantilever Beam. Comparison of Full and 
Reduced Integration Results for P = 40 B/in. 





Total Average 



Tip Deflection 

strain energy release rate 


Case 

w, in 

G lb. in/ in 2 

1 . 

Full Integration 

.0323 

1.129 

2. 

Reduced Integration 

0.06416 

4.780 


of t only (IREDU=1) 
yz 



3. 

Reduced Integration 
of all Shear Strains 

0.06595 

4.928 


( IREDU=3) 



4. 

Beam Theory 

0.06540 

4.915 
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along the crack-front. As should be expected the strain energy release rate 
has maximum value at the center and becomes negligibly small at the edge. 
Since the strain energy release rate decays to zero only in a small portion 
of the width at the edge, a finer mesh is necessary to get the accurate 
distribution of G. A finer mesh would result in a flatter curve in the 
central portion and a steeper decay at the edge. 

4.3 Double Cracked Lap Shear Specimen 
A symmetric double cracked lap shear specimen made of laminated com- 
posite is shown in Fig. 7a. A through-the-thickness delamination of length 
a=l" exists symmetrically on either side of the specimen. Since the geom- 
etry is symmetric, not much of bending deformation is expected. Also, 
because of the symmetry, only one (upper) half of the specimen needs to be 
analyzed. NAS3D was used to carry out the 3-D analysis of this specimen. 

The mesh in the Y-Z plane used for the FE analysis is shown in Fig. 7b. The 
3-D mesh is generated by repeating this mesh in every Y-Z plane for the 
subdivision in X-direction. Full integration was used. The results of the 
linear solution were compared with those obtained earlier by Raju (Ref. 9). 
The results of displacements and stresses were identical with Raju's for the 
identical meshes. The strain energy release rates along the crack-front 
obtained in these solutions are compared in Fig. 7b. Although the average 
values in the two solutions match some difference is observed in the actual 
distribution. This is due to the different techniques adopted to get energy 
per unit area from nodal forces and displacement. The present method of 

calculating consistent tractions and using them in G-calcul ation is expected 
to be more accurate than the earlier one. 
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Example of DCLS specimen, (a) configuration, and (b) FE mesh in y-z plane. The mesh division 
is repeated on each plane parallel to this plane to get 3-D mesh. 


4.4 Cracked Lap Shear (CLS) Specimen 

Cracked lap shear specimens are widely used in the study of mixed mode 
fracture related to debonds in adhesively bonded joints or delaminations in 
composites. One such specimen used for adhesive bond study is shown in Fig. 
9a. The adherends are of equal thickness and made of aluminum with Young's 
modulus of 10500 ksi and Poisson ratio of .33. The adhesive Young's modulus 
is 280 ksi and Poisson ratio is 0.4. This specimen is one of the specimens 
used in a recent ASTM round robin effort on the analysis of CLS specimen 
(Ref. 7). The analysis was conducted using NAS3D for the case of debond 
length of 1". Reduced integration was used. The linear solution for P=2500 
lb. was obtained. The distribution of strain energy release rates in 
various modes along the crackfront is shown in Fig. 9b. These results agree 
well with the results of an earlier 3-0 analysis by Lof for the ASTM round 
robin (Ref. 7). However, the total average strain energy release rate is 
much higher than those obtained in nonlinear and linear 2-D analyses. This 
is surprising; but, it is expected that the nonlinear analysis will give the 
strain energy release rate comparable to the 2-D solutions. Further, to 
check the 3-0 program plan strain conditions were simulated in the 3-D 
analysis by restraining the normal displacements of the edges of the speci- 
men. The results of this analysis were within 2 % of the 2-D plan strain 
linear analysis. 

During the course of nonlinear analysis of this specimen, it was found 
that the initial load and the incremental load values required to obtain 
convergence were very small. Thus, nonlinear solutions for any load values 
of practical significance would require large computation time. 

5. CONCLUDING REMARKS 

An algorithm for computation of stiffnesses required in geometric 
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front in the linear solution. 





Adhesive 



Example of CLS specimen, (b) configuration, (b) G-distributions along the crack front 





nonlinear analysis was developed and used to generate the computer program 
NAS3D. Selective reduced integration of shear strains was incorporated to 
make the element suitable for bending. Strain energy release rate calcula- 
tions were based on consistent tractions along the crack front. A provision 
for restarting the solution from an intermediate solution is also made. The 
code was tested by solving some simple problems having known solutions. 
Comparison of various computation times showed that savings by a factor of 
11-13 were obtained in vectorized computation of stiffness over the conven- 
tional scalar computation. Calculation of residuals was found to be another 
time consuming step and future efforts may be directed to make this more 
efficient. 

Successful use of any FE program for the structural analysis depends 
significantly on the ability of the analyst to predict qualitatively the 
response of the configuration. This insight, generally based on experience 
and simpler analyses, is particularly important for nonlinear analyses in 
which questions of convergence and uniqueness of solution and solution 
strategy have to be addressed. It is hoped that the present program coupled 
with some of this insight will prove to be a useful tool for the geometric 
nonlinear analysis. 
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APPENDIX 


INPUT DATA 

The required input data is described in this section. Where 
applicable, the maximum allowable values of the input parameters are noted. 

No. of 

Card Set Parameters Cards Format 

1. TITLE ( I ) » I = 1,60 3 20A4 

TITLE = TITLE OF PROBLEM 

2. OUTPUT, ANALYS, QUADRAT, ENERGY 1 5A8 

OUTPUT = Output option 

= XLONG for long output 

= SHORT for output (the nodal coordinates, element 
connectivity, and boundary conditions are not in the 
output) 

ANALYS = Type of analysis 

= XLINEAR for linear analysis 
= GNONLIN for geometrically nonlinear analysis 
= PNONLIN for materially nonlinear analysis 
= CNONLIN for combined geometric and material nonlinear 
analysi s 

QUADRAT = Integration option 

= REDUC for reduced integration 
= XFULL for full integration 

ENERGY = Option for strain-energy release rate calculations 
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Card Set 


Parameters 


No. of 
Cards 


Format 


= DO-G for G calculation 
= DONOJG for no G calculation 

3. IREDUO 1 1615 

IREDUO = Number of shear strains to be reduced integrated 

U for T yz , 2 for r yz and t zx , 3 for T yz , and r xy ) 

4. ITSTEP, NCYCLE, IMAX 1 315 

ITSTEP = Number of steps in the incremental loading minimum = 1, 
maximum = 30 

NCYCLE = Number of iterations between updates of stiffness 
matrix 

IMAX = Maximum number of iterations allowed before terminating 

5. ACCURACY 1 FI 0.3 

ACCURACY = Maximum residual allowed in converged solution as a 
fraction of applied load 

6. NN, NE, NRN, NLX 1 1615 

NN = Number of nodes in the FE model 

NE = Number of elements in the FE model 

NRN = Number of nodes with a restrained degree of freedom 

NLX = Number of nodes in the X-direction (direction normal to 

the 2-0 mesh) (minimum = 1, maximum = 10) 

7. X(I), I = l.NLX NLX/8 f 8F10.5 

X( ) = Coordinates in X direction for mesh division 

8. Nodal Coordinates 

x-coordinate 
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Card Set 


Parameters 


No. of 
Cards 


Format 


XX, N( I) = 1,13 * E10.4, 1315 

XX = coordinate 
XX = coordinate 

N( ) = list of nodes with coordinate XX 
♦Input until all x-coordinates are 
specified. End x-coordinate data 
with a blank card. 
y-coordinate 

XX, N( I ) , I = 1.13 * E10.4, 1315 

♦Similar to input of x-coordinates 
9. I, I N ( I ) , JN( I ) , KN( I ) , LN( I ) NE 515 

I, IN, JN,KN,LN = Element number, four node numbers for element I. 

Nodes must be specified in a counterclockwise 
direction. 

10. K, NRL (3*K-2) , NRL (3*K-1), NRL (3*K) NRN 415 

K = Node number 

NRL (3*K-2) , NRL (3*K-1), NRL (3*K) = Constraints in X, Y 

and Z directions, 
respectively, at node 
K. 0 indicates no 
constraint 
1 indicates 
constraint 

Note: Do not include degrees of freedom involved in 
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No. of 

Card Set Parameters Cards Format 

multipoint constraints. Do include degrees of 
freedom with specified displacements. 

11. IANC 1 1615 

IANC = Number of additional nodal constraints on planes or 
1 ines. 

12. XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX 3* 2F20.8 

XMIN, XMAX = Range of x-coordinate 
YMIN, YMAX = Range of y-coordinates 
ZMIN, ZMAX = Range of z-coordi nates 

13. ICONX, ICONY, ICONZ 1* 1615 

ICONX = Restraint in x-di recti on 
ICONY = Restraint in y-di recti on 
ICONZ = Restraint in z-direction 
(=1, restrained j = 0, free) 

SKIP 14-17 IF ENERGY = OONOJG 

14. INP 1 15 

INP = Number of node sets used in virtual crack extension 
calculation (maximum = 15) 

15. NEGCAL(I), I = 1, ( INP+1) (INP+D/IG 1 1615 

NEGCAL = Element numbers for elements in 2-D mesh contributing 
to the nodal forces required for virtual crack 
extension. (See example in sketch below. Element 
numbers are circled.) 

^ Round off to next higher integer. 
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No. of 

Card Set Parameters Cards Format 



IF INP = 3, 

NEGCAL (1 to 4) = 2, 3, 4, 5 

NFGCAL (1 to 3) = 14, 13, 12 

NDGCAL (1 to 6) = 15, 19, 16, 20, 17, 21 

16. NFGCAL(I), I = 1, INP INP/16 1 " 1615 

NGGCAL(I) = Node numbers for nodes along which virtual crack 
extension forces are calculated. 

List according to distance from crack tip, with the 
crack tip node as the first one. 

17. NDGCAL ( I ) , I = 1, (2*INP) 2*INP/16 t 1615 

NDGCAL(I) = Node numers for the nodes used to calculate cracking 
opening and sliding displacements 
Repeat card sets 13-16 for each material group. 

Maximum number of material groups = 10 
End last group with blank card. 

18. J, XMATER(J), ELASOPT 1 15, A8 

J = Material group number 
XMATER = Material type 

= ELASTIC for linear stress-strain curve 
= ELPLAST for el asti c-perfectly plastic stress-strain 
^Round off to next higher integer. 
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Card Set 


Parameters 


No. of 
Cards 


Format 


curve 

= BLINEAR for bilinear stress-strain curve 
= RAMOSGO for Ramberg -Osgood stress-strain curve 
ELASOPT = YUNGMQD, when elastic constants of tne material are 
i nput 

= FULDMAT, when [D] is input columnwise 
(Use ELASOPT = YUNGMOD with the present version) 

19. EPROP(I), I = 1, 10 2 8E10.3 

EPROP = Contains the material elastic constants in the 
following order 

^1* ^2 » ^3 » ^2 3* ^3 1 » ®32* v 2 3* v 31’ v 12’ ® where 

0 = Ply angle with respect to x-axis 

20. YIELDS, ET, RO, ANM 1 5E10.3 

YIELDS = Yield stress 

ET = Tangent modulus for yielded bilinear material 

RO, ANM = Parameters defining Ramberg-Osgood stress-strain 

relation, e = ! + (?-) ANM 
E RO 

(a) If XMATER = ELASTIC, input YIELDS = ET = RO = 1.0 x 
10 21 , ANM = 10 

(b) If XMATER = BLINEAR, input proper YIELDS and ET and 
set RO = ANM = 0.0 

(c) If XMATER = RAMOSGO, input proper YIELDS, RO and 
ANM and set ET = 0.0 
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Card Set 


Parameters 


No. of 
Cards 


Format 


21. NEL1, NEL2, NELINC * 315 

NEL1, NEL2, NELINC = Loop parameters used to define elements in 

material group 
NEL1 = First element 
NEL2 = Last element 
NELINC = Loop increment 
e.g., 1, 50, 20 defines elements 1, 21, and 
31 to be in material group 
♦Repeat until all elements in group are defined. 

End card set 16 by specifying NEL1 = NEL2 = NELINC = 0 

22. OELLOAD(I) = 1, ITSTEP ITSTEP/8 + 8F10.3 

DELLOAD(I) = Scale factor for proportional load vector for load 
step I. Always specify DELLOAD(I) = 1.0 

23. NLN, NCD, NED 1 315 

NLN = Number of nodes with applied loads 

NCD = Number of multipoint constraints, max = 15 

NED = Number of specified displacements, max = 30 

24. K, FX, FY, FZ NLN 15, 3F10.3 

K = Node number 

FX, FY, FZ = Loads in x, y, and z directions, respectively 

25. K, KDF, URD NED 215, F10.3 


K = Node number 

KDF = Displacement direction, specify 1 for x direction 

specify 2 for y direction 


t 


Round off to next higher integer. 
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No. of 

Card Set Parameters Cards Format 

specify 3 for z direction 
URD = Magnitude of displacement 
SKIP 21-24 IF NCD = 0 

26. NMPR(I), I = 1, NCD NCD/16 1615 

NMPR(I) = Number of degrees of freedom involved in the It h 
multipoint constraint, max = 20 

27. ( ( ICDN( I , J) , J=l, NMPR(D), I = 1, NCD) NCD sets 1615 

( ICON ( I , J ) = Jth degree of freedom involved in the Ith multipoint 
constrai nt 

28. NZKV . 1 15 

NZKV = Number of multipoint constraints for which there is 
an applied load 

29. NKV, ATOT 

NKV, ATOT: ATOT is the non- zero load associated with the NKV 

set of constrained nodes 

30. NPANS 1 15 

NPANS = Number of planes/lines with applied uniform normal 
stress 

SKIP 31 if NPANS = 0 
Repeat 31-33 NPANS times. 

31. NDIR, ASTR 1 15, F20.8 

NDIR = Direction of applied stress 
(1 = x, 2 = 4, 3 = z) 

ASTR = Value of applied stress 
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No. of 

Card Set Parameters Cards Format 

32. XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX 3 2F20.8 

XMIN, WAX = Range of x coordinate for the plane/line 
YMIN, YMAX = Range of y coordinate for the plane/line 
ZMIN, ZMAX = Range of z coordinate for the plane/line 

33. NI, NF, NPCHX, NPCHY, NPCH2 1 1615 

NI = Lowest node nunber in the plane/ line 
NF = Highest node number in the plane/line 
NPCHX = Pitch in x-direction of node nunbers 
NPCHY = Pitch in y-direction of node numbers 
NPCHZ = Pitch in z-direction for node numbers 

34. IREDO 1 15 

IREDO = Restart parameter 

(0 = regular job, 1 = restart job) 

Note: The information required for restarting a job from an 

intermediate load step is written on TAPE7. This file should 
be saved and used for restarting the job. Thus, in the first 
run of the job IREDO is given as 0 and file on TAPE7 is saved. 
In the next run (for restarting at the point where the first 
run was over), IREDO is given as 1 and an earlier solution is 
read from TAPE7. 
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